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A two-dimensional rapidly rotating Bose-Einstein condensate in a harmonic plus quartic trap is 
expected to have unusual vortex states that do not occur in a pure harmonic trap. At a critical 
rotation speed flh, a central hole appears in the condensate, and at some faster rotation speed 
fig, the system undergoes a transition to a giant vortex state with pure irrotational flow. Using a 
time- dependent variational analysis, we study the behavior of an annular condensate with a single 
concentric ring of vortices. The transition to a giant vortex state is investigated by comparing the 
energy of the two equilibrium states (the ring of vortices and the giant vortex) and also by studying 
the dynamical stability of small excitation modes of the ring of vortices. 

PACS numbers: 03.75.Kk, 67.40.Vs, 31.15.Pf 

I. INTRODUCTION 

Since the creation of a single vortex in Bose-Einstein condensates in harmonic traps 0, 0] , experimental techniques 
have rapidly created arrays of several hundred vortices [3, Lil L2 111 • In these experiments, the condensate rotates at 
angular velocities fl that typically approach the high-rotation limit lu± , where the centrifugal effect from the rotation 
would cancel the radial confinement of the trap • 

For fl — > u>± in a pure harmonic radial trap, the combined effective radial trap vanishes and the Thomas-Fermi 
(TF) radius of the condensate diverges. Moreover, as Ho pointed out, the single-particle Hamiltonian in a rapidly 
rotating harmonic trap is analogous to that of a charged particle in a uniform magnetic field, and all particles should 
condense into the lowest Landau level (LLL). This insight has stimulated several experimental 0,^3 and theoretical 
[HIHGl works. 

An additional radial trap potential, stronger than harmonic, opens new possibilities in the study of rapidly rotating 
Bose-Einstein condensates, because it can provide confinement even for fl > uj±. A quartic potential is a simple choice 
0, 0| that has recently been realized experimentally with a blue-detuned laser directed along the axial direction of 
an elongated condensate |l5l [16J . Not only does the combined potential overcome the limit on the angular velocity 
f2, but it also enriches the physics of rapidly rotating condensates with the possibility of the interesting new vortex 
configurations. 

In such a combined potential, no quantized vortices appear in the condensate at sufficiently slow rotations fl <C lo±. 
With increasing fl (< u>±), the condensate passes through a sequence of states with many vortices that eventually 
form a large vortex lattice. These states are essentially the same as those in a pure harmonic trap. The new vortex 
phases begin to appear as the condensate continues to expand radially for ft > uj±. Specifically it develops a central 
hole at a critical angular velocity flh that depends on the interaction strength and the trap potential. For larger 
values of fl > O^, the condensate forms an annulus, whose mean radius expands and whose width continues to shrink 
with increasing fl. Finally a "giant" vortex state with pure irrotational flow is expected to appear at fl g , when all 
the singly quantized vortices disappear from the annulus. 

Previous theoretical work has either relied on numerical methods 0, or dealt only with weak interactions and 
small anharmonicities 0,H(j- In contrast, we here use a time-dependent variational analysis to study the behavior 
of an annular condensate with a single ring of vortices in the annulus in the TF limit. As Q increases beyond Oft, the 
width of the annular condensate becomes increasingly narrow, and above some value of angular speed fl (smaller than 
fl g ) the annular condensate can support only a one-dimensional array of vortices arranged in a single ring |l7ll2lll22| . 
The transition to a giant vortex state is investigated by comparing the energy of these two equilibrium states (a state 
with a single ring of vortices and a giant vortex state with pure irrotational flow) and also by studying the dynamical 
stability of small-amplitude normal mode excitations of the vortex array. Section ^ summarizes the basic procedure 
of the Lagrangian formalism and the analytical study in the TF limit. In Sec. IIIII we introduce our basic approximate 
condensate density profile, which is simple, yet incorporates the discrete character of the quantized vortices. We 
analyze the equilibrium (Sec. lIIll) and dynamical (Sec. lIV|) behavior of the condensate using variational methods, and 
the transition angular velocity fl g is evaluated numerically for various values of the interaction parameter and the 
trap parameter. 
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II. BASIC PROCEDURE 



A trapped Bose-Einstein condensate at zero temperature can be described with a macroscopic order parameter (the 
condensate wave function) <3> = (iff). It satisfies the time-dependent, nonlinear Schrodinger equation (Gross-Pitaevskii 
equation), 

1% 2 V 2 Anh 2 a s lTl2 T 

where a s is the positive s-wave scattering length, characterizing the repulsive interaction between the condensate 
particles. With an additional quartic component, the external two-dimensional trap potential V ir can be written as 



VtI = -Mu;l(r 2 +X- T )=-hu x - 2 - + A- r ), (2) 



where d± = y^h/McoZ is the harmonic oscillator length, r is the radial coordinate in two dimensions, and A is a 
dimensionlcss parameter that determines the relative strength of the quartic component. 

If the external trap potential rotates at an angular velocity O = Qz, it is necessary to work in a co-rotating frame, 
and the time-dependent Schrodinger equation becomes 

^-^™ + ^i*r-*-n.L*, ,3, 

where L = r x p = — ihr x V is the angular momentum operator. 

Frequently, instead of this Schrodinger equation, it is preferable to use an equivalent Lagrangian formalism, 

C[*] = T [*]-?[*], (4) 

where 

(•*£-£•) 

is the time-dependent part of the Lagrangian functional and 
T [*] = £' [*] - (j,N 

Ti 2 , l9 , ,o 27T?i 2 a, 



2M , V*| 2 + V tI |*| 2 + |*r - • - M 



dV 

2M 1 1 " ' ' 1V1 ' ' ' 
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(6) 



is the free-energy functional, incorporating the constraint of fixed particle number N = J dV Y$>\ with the chemical 
potential [i as a Lagrange multiplier. 

For simplicity, we consider an effectively two-dimensional system, uniform in the z direction over a length Z . The 
condensate wave function can be chosen as * = a/ A? '/Zip, and the normalization condition becomes 

l = Jd 2 v\^\ 2 . (7) 

In dimensionless form, with d± and lu± as scales for length and frequency, the free-energy functional (per particle) 
becomes 



T = d 2 r 



\ I W| 2 + \ (r 2 + Ar 4 ) M 2 + | |V| 4 + ty*0 txV^-m M 2 



(8) 



where g = AnNa s /Z is a dimensionless coupling constant, characterizing the interaction strength between the con- 
densate particles. The physics of the condensate is now determined by two dimensionlcss constants: the quartic trap 
strength A and the interparticle strength g. 

Using the replacement ip = ^Jne lS ' , we can write the first term on the right hand side in Eq. JSJ) as |V?/;| 2 = 
IVv 7 "! + IVS^ 2 \ip\ 2 ■ In our units, the dimensionless particle velocity field is given by v — VS* and, in the Thomas- 
Fermi approximation, the curvature of the density V-y/n is neglected. In this limit, variation of the free energy T 
with respect to I yields the familiar TF density 

g\^)f = (i + Ct ■ r x v - - (v 2 + r 2 + Ar 4 ) . (9) 



Given the velocity v, we can determine /x with the normalization condition, Eq. Q. 
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A. Uniform vortex lattice with central hole 



In the limit of a dense uniform vortex lattice with dimcnsionless areal density Q/tt |23j, the total velocity field 
arising from each quantized vortex in the annular region i?< < r < R > can be approximated by an integral, yielding 

v a (r) - n (r - ^ \ 4>. (10) 

The first term is the anticipated solid-body rotation, but the second term reflects the presence of the central hole in 
the annulus. To better mimic solid-body flow that minimizes the free energy in the annulus, it is necessary to add 
the irrotational flow Vi rr (r) = (flR^/r) 0, arising from the phantom vortices in the hole. Then the total superfluid 
flow in the annular region becomes the sum of the contribution v a from vortices in the annulus and the irrotational 
flow Vi rr from the phantom vortices, v = v a + vi rr = $1 x r = Qr<fi, which is the expected solid-body flow. 
Substituting this velocity field, the TF density in Eq. © is simplified to 

ff |0| 2 =/i+i[(O 2 -l)r 2 -Ar 4 ]. (11) 
This density profile predicts the appearance of a central hole in the condensate at angular speed f£/j that is given by 

M 

^ = 1 + 2VA^V /3 . (12) 
For Q > £2 the squared TF inner and outer radii of the annulus satisfy the simple relations 

1/3 



B. Irrotational flow in annular region ("giant vortex") 

The superfluid velocity field of pure irrotational flow is simply v = (y jr) <p, where v is the quantum number of 
the central circulation [y is an integer.) At sufficiently fast angular velocity, v will be large, and we can treat v as a 
continuous variable. The TF density for a giant vortex is then given by 

g\4>\ 2 =fi-U(r 2 ), (14) 

where fx = [i + Slv and 

U (X) = \ + X + Xx2 ) (15) 

can be considered an effective potential that combines the centrifugal barrier and the original trap potential (here, 
x = r 2 ). 



III. VARIATIONAL ANALYSIS 



As explained in the Introduction, the width of the condensate becomes narrower with increasing angular velocity, 
and at some point, the annular condensate can support only a one-dimensional array of vortices arranged in a single 
ring. We use a variational analysis to study a condensate with irrotational flow around a central hole (the giant vortex 
state) and a condensate with a single ring of N a singly quantized vortices in the annulus (the N a ^ state). 
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A. Simple parabolic density approximation 



Although the TF density profile of a giant vortex state in Eq. I|14(l is simple and manageable, when we include 
vortices in the annulus, it becomes complicated and almost untractable for variational analysis. We need a simpler 
density approximation, yet it should contain the same physics as the full system. 

From the TF density profile of a giant vortex state <|14(l . we know that the condensate is confined in an annular 
region bounded by inner and outer radii, and its density profile is a function of r 2 . As a simple model of the density 
profile for an annular condensate, we use the following parabolic functional form. 

gn^g\^\ 2 ^A(X > -r 2 )(r 2 -X < ), (16) 

where X < and X > are variational inner and outer squared radii of the condensate, respectively, to be determined 
such that they minimize the free energy of the condensate. A positive constant A > is also introduced to ensure 
the normalization condition. 

The comparison to the TF density profile of a giant vortex state is shown in Fig. for A = 1/2 and g — 1000, from 
fi = 4 to fi = 9. The optimized density profile from Eq. (|l(j[) is plotted in solid lines and compared to the TF profile 
(dotted). For easy comparison, the difference of two profiles (the parabolic density and the TF density) is also shown 
in dashed lines. Even at low fi values (for example fi = 4), this simple approximation shows good agreement - when 
the difference is compared to the maximum condensate density, it is < 5%. The approximation gets better as the 
external rotation speed fi becomes faster, and at fi = 6, the two density profiles are almost the same. Generally, the 
parabolic density profile has smaller inner and outer radii, but a slightly higher peak. 

We start from this trial parabolic density profile, along with an appropriate velocity field, and calculate the free- 
energy functional T . Later we repeat the procedure with a different improved model. As in |22j, the irrotational flow 
in the annulus can be approximated by the central circulation v. The total velocity field given by the phase S is 



v = VS= -0 + ^v fc , (17) 



where Vfe is the contribution from the A;th vortex located at = (rk,4>k) as one of a single ring of vortices in the 
annular region, 

r-r k 

^k=zx- (18) 

|r-r fc | 

In terms of the trial wave function ip = y/ne , the free-energy functional J- is now characterized by a set of 
parameters, T ' = T ' (fi, u, X < , X > , A, {r^}). For a given fi, we choose a value of N a (since the number of vortices N a 
is not a continuous variable) and we then minimize T with respect to all parameters other than /i, 

^_dT_ _ dT_ _ dT_ _ dT_ _ dT_ _ dT_ 
~ dv ~ dX< ~ dX> ~ dA ~ dr k ~ d<p k ' 
along with the normalization condition that determines fi with Eq. Ifllty. 

1 = l^(X > -X < f. (20) 
o g 

Once we minimize T and have all the values of the other parameters for given (fi, N a ), then we can compare the 
energy £' = T + fj, among different N a vortex states to determine the lowest energy state (ground state) for a given 
fi. 

Table [I] shows (meta)stable states for some typical values of fi at A = 1/2 and g = 1000. Here we only consider 
N a 7^ states, omitting the giant vortex state from the table. R = (VA> + \/A<) /2 is the mean radius of the 
condensate, R r is the variational radius (r^) of a symmetrically spaced ring of vortices in the annular region, b is the 
intervortex spacing b = 2itR r /N a , and d = V / A > — \/X < is the width of the condensate. The conventional (usual) 
healing length £o = 1/ V^9 n o 1S determined from the maximum density of the condensate uq = [A/Ag) (X> — X<) . 
The healing length £o in Table [I] is nearly constant, since the maximum density no remains almost the same as seen 
in Fig. ^ For comparison, the result from a numerical analysis from Ref. is shown in Table ITU 

As seen in Table [ij our variational study with the simple parabolic approximation in Eq. (|lfc>|) predicts that the 
intervortex spacing b decreases with increasing fi, although somewhat slower than the dependence fi -1 / 2 associated 
with a uniform vorticity. In contrast, the numerical study of the full GP equation |22| found that the intervortex 
spacing b increases slowly with increasing external rotation fi. Evidently, the simple parabolic approximation tends to 
overestimate the number of vortices in equilibrium. This qualitative feature must be corrected if we expect a realistic 
result from our model. 
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FIG. 1: The simple parabolic density (E) and the TF density JUJ profiles for A = 1/2 and g = 1000, at = (4, 5, 6, 7, 8, 9). 
The parabolic profiles are shown in solid lines, the TF profile in dotted lines, and, for comparison, the difference between the 
two is dashed. In all graphs, the horizontal axes (radius r) are drawn in the same scale, and we can easily see that the mean 
condensate radius becomes bigger and the width of the condensate becomes narrower with increasing Q. 



Q 


N a 


R 


Rr 


b 


d 


£o 


5 


45 


4.82 


4.86 


0.678 


1.88 


0.138 


6 


57 


5.87 


5.89 


0.650 


1.51 


0.138 


7 


71 


6.90 


6.92 


0.612 


1.29 


0.139 



TABLE I: Result of a simple parabolic approximation for A = 1/2 and g = 1000: Here N a is the optimized vortex number in 
a ring for given Q,. The mean radius is R = U/X> + s/X<) /2, R r is the variational radius of a ring of vortices, b = 2nR r /N a 

is the intervortex spacing, d = \/X> — \/X< is the width of the condensate, and £o = l/V2<?^o = \f%/A (X > — X<) _1 is the 
healing length. 



B. Modified parabolic density approximation 

We can improve the simple parabolic model by including the contribution of the nonuniform density near the core 
of each vortex to the total free energy. This effect makes the addition of a vortex in the annular region more costly and 
thus tends to reduce the number of vortices in equilibrium, increasing the intervortex spacing b. This core contribution 
becomes more significant as the radial trap rotates faster and the number of vortices increases. 

As a model for the vortex core structure, we choose a simple linear model [TTl liH l2~i| . 



i> = Vne iS F, 



(21) 



6 



n 


N a 




_Rtf 


frmcas 


^mcas 


^TF 


5 


37 


4.83 


4.79 


0.821 


2.32 


2.06 


6 


44 


5.80 


5.86 


0.828 


1.98 


1.68 


7 


51 


6.85 


6.89 


0.844 


1.76 


1.43 



TABLE II: Result of a numerical analysis for A = 1/2 and g = 1000 22]: N a here is determined by solving the time-dependent 
GP equation. i? mca s, fcmcas and d mcas are the mean radius, the intervortex spacing and the width of the annulus from the 
numerical calculation, respectively. Note that the solution of the full time-dependent GP equation does not have well-defined 
inner/outer radii where the condensate density vanishes. Thus, the determination of d mcas depends on the numerical technique 
used to identify the radii, and only the qualitative comparison is meaningful. For comparison, the TF estimates from Eq. I13H 
are also shown: Rtf = (R> + R<) /2 and c£tf = R> — R<- 



F = l[fk, (22) 



k 

where n is the parabolic density in Eq. I|16(l and 

f =/ 1 ' if |r-r fc |>£ 

Jk \|r-r fc |/£, if|r-r fc |<£ 1 ' 

is a modified condensate density due to the presence of the vortex. Here we assume that the cores of the vortices do 
not overlap. We use the symbol F for notational convenience to denote the product. Note that £ is now a variational 
parameter that determines the size of the vortex core; it should be distinguished from the previous healing length 
£o = l/V^S^o that is determined directly from the maximum condensate density uq. 

With this modified density profile, the new normalization condition becomes, up to corrections of order £ 2 , 

i«~(*>-*<) a -^£»M). (24) 

y k 

where £ 2 is assumed small (measured in units of d\). We can recalculate the free-energy functional T , using 

VV> - {iVSF^i + VF^ + FVV^) e lS , (25) 



|W| 2 = \VS\ 2 F 2 n 



\VF\ 



F 2 I VVnl + ZFy/nVF ■ Vy/ri. 



(26) 



If we again neglect the terms involving V-y/n, we find 
./' = I d 2 r 

k J 





- {v 2 + r 2 + Ar 4 ) n — (Qrv^ + fi) n + ^ 



21 



- (v 2 +r 2 + Ar 4 ) - (f>r« + pL) 



n (\f k \ -1 



(27) 



+ \Jd 2 v{\VF\ 2 

The first term has the same functional form as the free energy of the simple parabolic density approximation in the 
previous subsection; in contrast, the last three terms come from the core modification in the density. The contribution 
to the free-energy integral from the last term in Eq. Ij26(l approximately averages to zero, since each V/& is cylindrically 
symmetric with respect to the vortex position and also because V^n f=a near the vortex position. 

As in the simple parabolic density approximation, we can again find the lowest energy (ground) state. The only 
difference is that we have one more variational parameter £, determined by the condition = dJ-/d^. 
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FIG. 2: Energy (per particle) £' of various N a states for A = 1/2 and g = 1000, at = (4, 5, 6, 7, 8, 9). Note that, for Q = 9, 
we cannot find consistent solutions for 1 < N a < 18. 



Figure |3 shows the energy £' in the rotating frame at several distinct values of Q = (4,5,6,7,8,9) of various N a 
vortex states that minimize T with respect to the other parameters. At Q — 4, for example, a giant vortex state 
(N a = 0) is clearly unstable and N a — 26 state is the lowest state among the states with a single ring of vortices. A 
faster external rotation (£1 = 5 or 6) makes a giant vortex state a local minimum and therefore metastable although 
the true minimum has N a 0. In the Q = 7 figure, a giant vortex state is the ground (lower-energy) state and the 
N a = 45 state is now only a local minimum and therefore metastable. Finally in the Q = 9 figure, only a giant vortex 
state is stable, for the local minimum at finite N a no longer appears. 

The ground-state (or metastable-state) energy values that we find through this variational procedure for 4 < f2 < 7 
with the modified parabolic density approximation agree well with the numerical result from the full GP equation 
[22|. We find that in general our model has lower energy (more negative) with a small discrepancy (< 5%), but 
it is probably expected because our TF model ignores the (positive) energy contribution of the condensate density 
curvature l|26l) . In this sense, it may be more meaningful to compare the energy differences between the irrotational 
giant vortex state (N a — 0) and the annular N a ^ state (or the vortex lattice state in [l^), through which the TF 



f 



increases at an almost constant rate 



effect effects should somewhat cancel out. In our model, AS' = ^ Na ^ — ^ Na 
of ~ 0.73 per Q, starting from a negative value; in contrast, although not as constant as in our model, AS' in the 
numerical result (2^ increases at a slower rate of ~ 0.64 per £1. 
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n 


N a 


R 


Rr 


b 


d 


c 


£o 


5 


33 


4.82 


4.80 


0.915 


1.83 


0.248 


0.130 


6 


40 


5.87 


5.85 


0.918 


1.53 


0.236 


0.131 


7 


45 


6.90 


6.87 


0.959 


1.31 


0.226 


0.131 



TABLE III: Result of the modified parabolic density approximation for A = 1/2 and g = 1000: £ is the variational vortex 
core size and the other parameters are defined as in Table The conventional healing length £o is slightly different from 
Table because of the somewhat different parameters. While the optimized vortex number N a is overestimated in the simple 
parabolic model (Table HJ compared to the result of the numerical analysis (Table HTj, the modified parabolic model instead 
underestimates N a by 4~6. Importantly, the spacing between the adjacent vortices b increases slowly with increasing Q, in 
qualitative agreement with the numerical finding from |22j. 



Table IIIII shows results from the modified parabolic density approximation for some typical values of tt. In this 
modified model, the energy of one vortex includes the core contribution. As we have anticipated, the core modification 
in the density reduces the number of vortices in equilibrium relative to the simple parabolic approximation. In fact, 
the values are somewhat underestimated by 4~6 when compared to the full numerical solution of 22], and the 
intervortex spacing b now increases with increasing ft, in agreement with the numerical result |22|. The variational 
core size £ is also shown in the table. While the conventional healing length £o is basically determined by the balance 
between the kinetic energy term and the interatomic interaction term in the free energy, the variational core size £ is 
determined by the velocity field (the phase contribution to the kinetic energy), the confining trap geometry and the 
interaction terms within the TF limit. When compared to the conventional healing length £rjj we find that the ratio 
£/£o * (1-73 ~ 1.91). For reference, the experiments with a small quartic trap potential (A ~ 10~ 3 ) and no central 
hole HE3 found that £/£ « 1-94. 

Figure |31 shows characteristic external rotation frequencies for three anharmonic parameters A = (1/8, 1/4, 1/2) and 
a given interaction parameter g = (250, 500, 1000). Solid squares denote the external rotation frequencies below which 
only states with a finite number of vortices are stable (ground states) and giant vortex states are not even metastable. 
Between solid squares and circles, finite- N a vortex states are ground states and giant vortex states are metastable. As 
the external rotation becomes faster, at rotation frequencies given by solid circles, a giant vortex state has the same 
energy as the optimum finite- N a vortex state. A giant vortex state becomes the ground state for an external rotation 
frequency faster than that denoted by a solid circle. Finally at faster than the solid triangle, only the giant vortex 
state is stable and presumably the ground state (at least within this model). 

Figure0]shows the number of vortices in the optimum symmetric ring in the annular region of a stable or metastable 
condensate as a function of the external rotation f2. Interestingly, in most cases above some relatively large value of 
fl (although states with finite N a are no longer ground states), the number of vortices in equilibrium indeed decreases 
with the increasing external rotation until the N a states become unstable. 



IV. TKACHENKO-LIKE OSCILLATION MODES AND THE ASSOCIATED STABILITY 



By assumption, our equilibrium one-dimensional arrays of vortices are symmetrically spaced around a circle. From 
one of these equilibrium ground states with N a ^ 0, we can consider the small-amplitude normal mode excitations of 
the vortex array by considering small variations about the equilibrium positions of the vortices in the annulus. 

The time-dependent part T of the Lagrangian (C = T — T) can be calculated from 



T 



1A 

e 



rl-X < ) 2 (3X > 



2rl- X< )-^rl (X, 
z 9 



X< - 2ri 



(28) 



2 v kl rf + r| - 2rjr k cas(j> jk 



As in our previous study of the dynamics of a single ring of vortices in trapped Bose-Einstein condensates [21j 
introduce small variations a k and j3 k for radial and angular coordinates of the fcth vortex. 



we 



Oik, 



= F k + 0k 



(29) 
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FIG. 3: Characteristic external rotation frequencies for A = (1/8, 1/4, 1/2) and g = (250,500, 1000). For f2 below the solid 
squares only states with iV a 7^ are stable. Between solid squares and solid circles, giant vortex states are metastable and 
N a / states are ground states. Between solid circles and solid triangles, the ground state is now the giant vortex state and 
N a =fc states are only metastable. Above the solid triangles, only giant vortex states are stable. 



where (ry, is the equilibrium position of the kth vortex and <fP k — 2irk/N a . 
The Euler-Lagrange equations of motion in terms of these small variations are 



W^U_W«*U + W** )A , (30, 

V V d( t>kdr q ) V V drkdr i ) TV d( t )kdr '~ ' 



and 



(d 2 T \ ( d 2 T d 2 T \ ■ ( d 2 T \ 



d 2 T 



0k, (31) 



where 1 < q < N a . In equilibrium, the N a vortices are symmetrically placed, and the condensate has iV -fold 
rotational symmetry. As a result, a Fourier decomposition of a and (3 



N a -1 



will decouple the original problem into N a sets of 2 x 2 matrix problems, 

d 



d P s = M s P s , (33) 
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FIG. 4: Equilibrium number of vortices in a single ring in the annular region of a stable or metastable condensate as a function 
of the external rotation. 



where P s = ya s ,(3 s j. The Fourier index s is analogous to the wave number in a linear array. The matrix elements 

of M s are determined from Eqs. (f 3()|) and (|3 1 1> - 

By solving these matrix problems, we can study the dynamical stability of the vortex system. For a given equilibrium 
state, we determine the excitation spectrum as a function of the external rotation frequency Q. If a particular 
excitation frequency is real, the corresponding small amplitudes (a,f3) oscillate around zero, and the state is stable. 
If any excitation frequency is complex with a positive imaginary part, the corresponding mode will grow with time, 
and the state is dynamically unstable. 

We find that all frequencies of these oscillation modes are real for all annular states with N a ^ in equilibrium; 
thus the single ring of vortices in our annular condensate is dynamically stable under small oscillations of the vortex 
positions. Note that this conclusion requires a detailed dynamical analysis. In other cases, for example a single ring of 
vortices in a circular two-dimensional condensate, the ring becomes dynamically unstable for sufficiently large vortex 
number [2l| . 

Our approximate variational Lagrangian trial function treats only the vortex positions as dynamical variables. In 
particular, we do not allow overall density variations. Thus this analysis differs from that in p5|. where linearized 
rotational hydrodynamics is used to analyze the low- lying normal modes of annular condensates. Comparison of these 
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two approaches merits further study. 

V. TRANSITION TO A GIANT VORTEX STATE 

According to the above study, the transition from an annular state with N a ^ to a giant vortex state is determined 
from the energetic consideration of the equilibrium states. It is not determined from the dynamical stability of the 
small-amplitude normal mode excitations of the vortex lattice. 

In the equilibrium analysis, we can identify two characteristic rotational frequencies that can be related to the 
transition; one is the frequency at which the ground state changes from the annular N a ^ state to the giant vortex 
state (denoted by solid circles in Fig. OJ), and the other is the frequency above which the annular vortex state is 
unstable (denoted by solid triangles in Fig.OJ). Between these two frequencies, the N a ^ state and the giant vortex 
state are energetically stable or metastable, and dynamically stable. As Jackson and Kavoulakis |2£j pointed out, the 
existence of multiple (meta)stable states indicates possible hysteresis effect in the phase transition as the rotational 
angular frequency varies. But strictly, in the presence of some disturbing process/noise, small enough not to destroy 
the system yet large enough to overcome the energy barrier between the two states, the condensate seeks the global 
energy minimum state. Thus we choose to identify the transition frequency f2 s to the giant vortex as that where the 
change of ground states occurs (solid circles in Fig. - 

VI. DISCUSSION 

We have used a variational analysis to study a rapidly rotating, two-dimensional annular Bose-Einstein condensate 
in a harmonic plus quartic radial trap, concentrating on the states with a single ring of vortices in the annulus. We 
have used a modified parabolic density profile approximation that incorporates the density change near the core of 
each vortex. Several parameters, inner squared radius A<, outer squared radius A> and vortex core size £ are varied 
to minimize the free energy J 7 , leading to the stable state. 

The time-dependent variational analysis indicates that a transition to a pure irrotational giant vortex state is 
determined solely from energetic considerations, and cannot be associated with the dynamical stability of small 
oscillation modes of vortex positions in this model. We identify the transition frequency £l g to the giant vortex state 
by comparing the energy of two equilibrium states, taking £l g as the frequency where the ground (global energy 
minimum) state changes from N a ^ to N a = 0. 

We assume the validity of the TF approximation, ignoring the contribution of the condensate density curvature to 
the free-energy functional. Since the usual healing length £ and the variational vortex core size £ are small compared 
to the intervortex spacing b and the width of the condensate annulus d in all ranges of Q used in this study, it should 
be safe to use the TF approximation. Our analysis also argues against the different scenario for the transition to a 
giant vortex state at larg e Q, where d becomes too narrow (d ~ £) to accomodate a vortex in the annulus and the TF 
assumption breaks down |22j |. 

We study systems of various values of g = (250, 500, 1000) and A = (1/8, 1/4, 1/2). For the smallest value of g = 125 
at A = 1/2, we cannot find even a metastable state with a single ring of vortices; only a giant vortex state is stable. 
In contrast, the numerical study at the same value of g and A indeed found a transition from vortex arrays to a 
giant vortex. The difference may reflect our reliance on the TF approximation, which holds for larger g. For g = 1000 
at A = 1/8, our transition angular frequency to the giant vortex state is Q g = 4.36. A calculation based on a uniform 
vortex lattice for similar values of g = 3607T and A = 1/8 finds a considerably smaller value of the transition 
frequency Q g « 2.40. For g = 1000 at A = 1/2, we find the transition angular frequency fl g = 6.76 to the giant 
vortex state, while the numerical simulation |22J finds no transition even at f2 = 7. Presumably this discrepancy arises 
from our approximate variational analysis, whereas the full numerical study |22j relies on the full time-dependent GP 
equation. It may also reflect the finite grid size used in [2^1. In addition, the vortices enter from the outer edge in 
the imaginary-time formalism, and they might become trapped in a metastable state with a ring of vortices, instead 
of reaching a giant vortex with irrotational flow pr| . 

It would also be interesting to extend the present analysis to the case of a negative harmonic term, where the 
condensate is annular even for fl = j2(| [271 • 
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